home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_sequence.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.3 KB  |  155 lines
  1. /* specfunc/bessel_sequence.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_bessel.h>
  26.  
  27.  
  28. #define DYDX_p(p,u,x) (-(p)/(x) + (((nu)*(nu))/((x)*(x))-1.0)*(u))
  29. #define DYDX_u(p,u,x) (p)
  30.  
  31. static int
  32. rk_step(double nu, double x, double dx, double * Jp, double * J)
  33. {
  34.   double p_0 = *Jp;
  35.   double u_0 = *J;
  36.  
  37.   double p_1 = dx * DYDX_p(p_0, u_0, x);
  38.   double u_1 = dx * DYDX_u(p_0, u_0, x);
  39.  
  40.   double p_2 = dx * DYDX_p(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
  41.   double u_2 = dx * DYDX_u(p_0 + 0.5*p_1, u_0 + 0.5*u_1, x + 0.5*dx);
  42.  
  43.   double p_3 = dx * DYDX_p(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
  44.   double u_3 = dx * DYDX_u(p_0 + 0.5*p_2, u_0 + 0.5*u_2, x + 0.5*dx);
  45.  
  46.   double p_4 = dx * DYDX_p(p_0 + p_3, u_0 + u_3, x + dx);
  47.   double u_4 = dx * DYDX_u(p_0 + p_3, u_0 + u_3, x + dx);
  48.  
  49.   *Jp = p_0 + p_1/6.0 + p_2/3.0 + p_3/3.0 + p_4/6.0;
  50.   *J  = u_0 + u_1/6.0 + u_2/3.0 + u_3/3.0 + u_4/6.0;
  51.  
  52.   return GSL_SUCCESS;
  53. }
  54.  
  55.  
  56. int
  57. gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v)
  58. {
  59.   /* CHECK_POINTER(v) */
  60.  
  61.   if(nu < 0.0) {
  62.     GSL_ERROR ("domain error", GSL_EDOM);
  63.   }
  64.   else if(size == 0) {
  65.     GSL_ERROR ("error", GSL_EINVAL);
  66.   }
  67.   else {
  68.     const gsl_prec_t goal   = GSL_MODE_PREC(mode);
  69.     const double dx_array[] = { 0.001, 0.03, 0.1 }; /* double, single, approx */
  70.     const double dx_nominal = dx_array[goal];
  71.  
  72.     const int cnu = (int) ceil(nu);
  73.     const double nu13 = pow(nu,1.0/3.0);
  74.     const double smalls[] = { 0.01, 0.02, 0.4, 0.7, 1.3, 2.0, 2.5, 3.2, 3.5, 4.5, 6.0 };
  75.     const double x_small = ( nu >= 10.0 ? nu - nu13 : smalls[cnu] );
  76.  
  77.     gsl_sf_result J0, J1;
  78.     double Jp, J;
  79.     double x;
  80.     int i = 0;
  81.  
  82.     /* Calculate the first point. */
  83.     x = v[0];
  84.     gsl_sf_bessel_Jnu_e(nu, x, &J0);
  85.     v[0] = J0.val;
  86.     ++i;
  87.  
  88.     /* Step over the idiot case where the
  89.      * first point was actually zero.
  90.      */
  91.     if(x == 0.0) {
  92.       if(v[1] <= x) {
  93.         /* Strict ordering failure. */
  94.         GSL_ERROR ("error", GSL_EFAILED);
  95.       }
  96.       x = v[1];
  97.       gsl_sf_bessel_Jnu_e(nu, x, &J0);
  98.       v[1] = J0.val;
  99.       ++i;
  100.     }
  101.  
  102.     /* Calculate directly as long as the argument
  103.      * is small. This is necessary because the
  104.      * integration is not very good there.
  105.      */
  106.     while(v[i] < x_small && i < size) {
  107.       if(v[i] <= x) {
  108.         /* Strict ordering failure. */
  109.     GSL_ERROR ("error", GSL_EFAILED);
  110.       }
  111.       x = v[i];
  112.       gsl_sf_bessel_Jnu_e(nu, x, &J0);
  113.       v[i] = J0.val;
  114.       ++i;
  115.     }
  116.  
  117.     /* At this point we are ready to integrate.
  118.      * The value of x is the last calculated
  119.      * point, which has the value J0; v[i] is
  120.      * the next point we need to calculate. We
  121.      * calculate nu+1 at x as well to get
  122.      * the derivative, then we go forward.
  123.      */
  124.     gsl_sf_bessel_Jnu_e(nu+1.0, x, &J1);
  125.     J  = J0.val;
  126.     Jp = -J1.val + nu/x * J0.val;
  127.  
  128.     while(i < size) {
  129.       const double dv = v[i] - x;
  130.       const int Nd    = (int) ceil(dv/dx_nominal);
  131.       const double dx = dv / Nd;
  132.       double xj;
  133.       int j;
  134.  
  135.       if(v[i] <= x) {
  136.         /* Strict ordering failure. */
  137.     GSL_ERROR ("error", GSL_EFAILED);
  138.       }
  139.  
  140.       /* Integrate over interval up to next sample point.
  141.        */
  142.       for(j=0, xj=x; j<Nd; j++, xj += dx) {
  143.         rk_step(nu, xj, dx, &Jp, &J);
  144.       }
  145.  
  146.       /* Go to next interval. */
  147.       x = v[i];
  148.       v[i] = J;
  149.       ++i;
  150.     }
  151.  
  152.     return GSL_SUCCESS;
  153.   }
  154. }
  155.